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ABSTRACT 

The collapse of dense cores with metallicities — IZq is studied by hydrodynamical calculations 
coupled with detailed chemical and radiative processes. For this purpose, we construct a simple 
chemical network with non-equilibrium reactions among 15 chemical species, H + , e, H, H2, D + , D, 
HD, C + , C, CO, CO2, O, OH, H 2 0, and O2, which reproduces the abundance of important molecular 
coolants by more detailed network very well. Starting from the initial density of 10 4 cm~ 3 , the evolution 
is followed until the formation of a hydrostatic protostar at the center ~ 10 21 cm -3 . In a lower- 
metallicity gas cloud, the temperature during the collapse remains high owing to less efficient cooling. 
After the cloud core becomes optically thick to the thermal emission by dust, temperature evolution 
at the center converges to a single trajectory, except for cases with metallicity < 10~ 6 Z Q , where the 
temperature remains slightly higher than in higher metallicity cases even after becoming optically 
thick to thermal radiation by the H2 collision-induced emission. The protostellar masses at their 
formation are a few 1O _3 M0, being slightly higher for cases with < 10~ 6 Zq. Using the temperature 
evolution at the center as a function the density, we discuss the possibility of fragmentation during 
the dust-cooling phase. The critical metallicity for the fragmentation is 10~ 5 Z Q assuming moderate 
elongation of the cloud cores at the onset of this phase. From the density and velocity distributions 
at the time of protostar formation, we evaluate the mass accretion rate in the subsequent accretion 
phase. The accretion rate is larger than the Shu accretion rate for the inside-out collapse from an 
initially static cloud ~ c 3 /G, where c s is the sound speed in the prestellar gas, by about a factor of 
ten owing to more dynamical nature of the collapse. Using these accretion rates, we also calculate 
the evolution of the protostars under the assumption of stationary accretion flow. For > 10~ 4 Zq, 
we succeed in following their evolution well after the arrival to the main-sequence phase. For lower 
metallicity cases, however, owing to too high accretion rates > a few 10 -3 -M©/yr, the total luminosity, 
which consists of contribution from accretion and internal one, reaches the Eddington limit, thereby 
rendering the stationary accretion impossible for > lOOAf©. Finally, we discuss possible suppression of 
fragmentation by heating of the ambient gas by protostellar radiation, which is considered important 
in the contemporary star formation. We argue that it is negligible for < 1CU 2 Zq. 
Subject headings: stars: formation, stars: Population II 



1. INTRODUCTION 

The first stars in the universe are thought to be formed 
from a chemically pristine gas. Recent theoretical studies 
and cosmological simulations suggest that these stars are 
very massive, perhaps even greater than 100 solar masses 
(Bromm, Coppi & Larson 1999; Abel, Bryan & Norman 
2002; Yoshida, Omukai & Hernquist 2008; Ohkubo et 
al. 2009). Stars in the solar neighborhood are, on the 
other hand, typically of low-mass with the initial mass 
function peaking around ~ 0.3M Q (e.g., Kroupa 2002; 
Bastian et al. 2010). Although the characteristic mass of 
old halo stars - Population II stars- is not well known, a 
substantial fraction of them were formed with low masses 
and thus are still observed in the Galaxy. Clearly, there 
must have been a transition in the characteristic stellar 
mass from massive to low-mass ones. 

Stars are formed under rather different conditions in 
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the early universe and in the local universe. For example, 
metallicity of the interstellar gas, strength of a radiation 
field, existence of cosmic rays and a magnetic field, are all 
thought to be critical quantities that affect the process 
of star formation. Gas metallicity appears to be among 
the most important quantities that affects the thermal 
evolution of a prestellar gas and hence might determine 
the characteristic stellar (at least gas cloud) mass (e.g., 
Bromm et al. 1999; Omukai 2000; Schneider et al. 2002). 

Observationally, no clear differences are found in the 
stellar initial mass function (IMF) for regions with dif- 
ferent metallicities in the local universe (e.g., Kroupa 
2002; Bastien, Covey, & Meyer 2010). However, the dif- 
ference of metallicity in the examined regions is only a 
factor of a few, and thus it remains unclear if the char- 
acteristic stellar mass can change drastically with much 
lower metallicity. Interestingly, there is an observational 
suggestion that very metal-poor stars in the halo had 
larger characteristic mass (Lucatello et al. 2005). The 
carbon-enhanced metal-poor stars are considered to be 
secondary, i.e., lower-mass, stars in binaries. The other 
star in the binary is typically an intermediate-mass star, 
whose mass is a few solar masses. The fraction of carbon- 
enhanced stars increases with decreasing metallicity, in- 
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dicating a higher fraction of intermediate-mass stars than 
in the solar neighborhood. 

Observations show that, in the Galaxy, the shape of 
dense-core mass-function resembles that of the stellar 
IMF (e.g., Mottc ct al. 1998). This fact indicates 
that the characteristic stellar mass is originally imprinted 
through the formation process of the dense cores, e. g., 
fragmentation of star-forming gas. According to recent 
detailed numerical simulations, the frequency of fragmen- 
tation depends strongly on the temperature evolution of 
gas. Fragmentation occurs frequently as long as the effec- 
tive ratio of specific heat 7 = dlogp/dlogp, where p is the 
pressure and p is the density of the gas, is less than unity, 
i.e., the temperature decreases with increasing density 
during the collapse. With 7 > 1, on the other hand, 
the fragmentation hardly occurs (Li ct al. 2003). The 
moment when 7 exceeds unity, i. e., at a temperature 
minimum, marks the end of the fragmentation epoch. 
The Jeans mass at this instant is imprinted as a char- 
acteristic fragmentation scale. Those numerical results 
can be interpreted by dynamical fragmentation theory 
of non-isothermal gas (Larson 1985), which predicts that 
the core mass-scale is set by fragmentation of filamentary 
clouds into round cores, which occurs when the isother- 
mality breaks. In the first star formation, dense cores 
are indeed produced around the temperature minimum 
at ~ 10 4 cm~ 3 , where the H 2 -cooling rate saturates owing 
to the collisional de-excitation. Also, for the present-day 
star formation, the Jeans mass at the dust-gas coupling, 
where the temperature becomes minimum, falls on the 
correct range for the characteristic stellar mass (Larson 
1985; 2005, Whitworth et al. 1998). 

Prestellar evolution of both present-day and primor- 
dial clouds has already been modelled in detail. In a 
primordial gas cloud, important processes are the chem- 
istry and radiative processes of H2, which are simple 
enough to incorporate into numerical simulations in a di- 
rect, self-consistent manner (Omukai & Nishi 1998; Abel, 
Bryan & Norman 2002; Yoshida, Omukai & Hernquist 
2008). The situation is clearly different for the present- 
day star formation. Quite a large number of radiative 
and chemical processes are operating in a metal-enriched 
gas cloud, e.g., cooling by atomic and molecular line 
emission, heating by the dissipation of turbulence and 
by external radiation fields, in particular, in a diffuse 
medium. At first glance, this complexity makes it ex- 
tremely hard to model the presente-day star formation 
process. In addition, once a central protostar turned on, 
the heating of ambient gas by its radiation becomes im- 
portant (Krumholz 2006; Bate 2009). However, consid- 
ering the fact that the cooling by dust grains dominates 
other radiative cooling processes after the formation of 
dense cores (> 10 4 cm~ 3 ) and before the protostar forma- 
tion, a simplification is possible: the prestellar collapse 
of dense cores can be modelled by considering only the 
radiative transfer of dust thermal emission. This makes 
the chemistry and radiative process effectively decouple 
(Gaustad 1963; Hattori, Nakano, & Hayashi 1969), and 
thus detailed calculations have been done since long ago 
(Larson 1969; Narita, Nakano, & Hayashi 1970; Winkler 
& Newman 1980; Masunaga & Inutsuka 2000). 

In Pop II star formation, the main difficulty lies in the 
coupling of chemical reactions and radiative cooling pro- 
cesses: one needs to solve chemical reactions to know 



the abundances of coolant species accurately. Miki & 
Nakano (1975) and Yoshii & Sabano (1980) studied the 
thermal evolution of collapsing low-metallicity clouds us- 
ing one-zone models. Similarly, Bromm & Loeb (2003) 
and Santoro & Shull (2006) examined the role of fine- 
structure line cooling in a low-metallicity cloud. These 
studies, however, do not follow chemical reactions but as- 
sume that metals in the gas phase remain in the atomic 
phase. Omukai (2000) and Omukai et al. (2005; here- 
after O05) solved the detailed chemical reactions as well 
as radiative processes. They identified the importance of 
cooling by CO, OH and H 2 6 for metallicity Z < 1O _1 Z . 
There still remains overall uncertainties, however. The 
dynamics is affected by the thermal evolution via pres- 
sure force, so the treatment in the one-zone models has 
not yet been justified. Three-dimensional hydrodynam- 
ical calculations were also performed by several authors 
(Bromm et al. 1999; Smith & Sigurdsson 2007; Smith 
et al. 2009). A slightly different approach is to use a 
prescribed effective cquation-of-statc in hydrodynamical 
calculations. In such studies, temperature evolution is 
given solely as a function of density, based on the re- 
sult of one-zone models (Clark, Glover & Klessen 2008; 
Tsuribe & Omukai 2006, 2008, Machida et al. 2008; 
2009). Jappsen et al. (2007; 2009a,b) included chemical 
reactions self-consistently in their hydrodynamical cal- 
culations, but followed only early evolution phase with 
density < 10 4 cm -3 . None of these calculations provide 
a complete picture of prestellar gas collapse with hydro- 
dynamics, chemistry, and radiative transfer. 

In this paper, we study the collapse of star-forming 
dense cores with different metallicities using radiation 
hydrodynamics calculations including chemical reactions 
for important molecular coolants. We follow their evolu- 
tion up to the formation of protostars at the center under 
the assumption of spherical symmetry. The structure of 
the cores at the moment of the protostar formation sets 
the initial conditions for the subsequent accretion phase 
of the protostar. In particular, it determines the proto- 
stellar accretion rate. Using the obtained accretion rates, 
we then calculate the evolution of protostars until they 
start hydrogen burning and arrive at the main sequence 
phase. 

This paper is organized as follows. In Sec. 2, the nu- 
merical method is described. We present the prestellar 
evolution of dense cores in Sec. 3 and then the proto- 
stcllar evolution in the accretion phase in Sec. 4. Sec. 
5 is for discussion. Finally, we summarize the paper in 
Sec. 6. In Appendix A the chemical network solved is 
presented. Appendix B contains the fitting parameters 
for the OH and CO cooling functions. 

2. METHOD OF CALCULATION 

We calculate the gravitational collapse of the dense 
cloud cores with metallicities Z = — 1 Zq by solving 
the radiation transfer and chemical reactions as well as 
hydrodynamics under spherical symmetry. Starting from 
the initial density ~ 10 4 cm~ 3 at the center, we follow the 
evolution until shortly after the formation of a protostar 
at - 10 21 cm- 3 . 

2.1. Hydrodynamics 

We solve hydrodynamics using the Lagrangian differ- 
ence method with the standard von Neumann- Richtmyer 
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viscosity (see, e.g., Mihalas & Weibel-Mihalas 1984). 
The results presented here are obtained with 200 grid 
points logarithmically spaced at the initial moment, with 
the innermost cell containing a mass of 1O _6 M0. The 
outer boundary is fixed at its initial position. To see the 
numerical convergence of our results, we also carry out 
calculations for the same setting but with the number 
of the cells halved. We have confirmed that the results 
agree well. 

2.2. Cooling Processes 

We include the same cooling processes as in O05 with 
some updates of the coefficients. The total cooling rate 
A consists of three processes, i.e., cooling via atomic and 
molecular line photons, cooling by continuum radiation, 
and the chemical cooling / heating: 

A = An ne + A cont + A c hem- (1) 

When the chemical heating exceeds the cooling, the last 
term A c h em becomes negative. 

2.2.1. Line-emission Cooling 

The line-cooling processes included are those by fine- 
structure transitions of CII, CI, and 01, and rotational 
transitions of H 2 , HD, CO, OH, and H 2 0. For H 2 , vi- 
brational transitions are also considered. Atomic and 
molecular line photons are trapped in the core when the 
optical depth is large. This effect is taken into account by 
reducing the cooling rate by a emission line multiplying 
the escape probability 

1 - e ~ TL 

/3l = —e-^, (2) 

where tl is the optical depth of the line and tq is that 
of continuum. The optical depth is evaluated with the 
radially-integrated column density 

N x (r) = f y{X)nndr', (3) 

J r 

where y(X) is the concentration of the emitting species 
X, and nn is the number density of hydrogen nuclei. 
The above integral is over the radial interval where the 
velocity shift is within the line width; \v(r') — v(r)\ < 
(v D (r r ) + v D (r))/2. In Omukai & Nishi (1998), transfer 
of H2 lines is solved not only by using correct averaging 
over solid angles but also by resolving frequency distri- 
butions within lines. By comparing the results with the 
same initial condition as Omukai & Nishi (1998) but with 
our approximations for line radiation transfer, we have 
confirmed that the results agree very well. 

The coefficient for CO, OH and H 2 line transitions 
are updated from O05. For H 2 0, we use the cooling 
rate given by Neufeld & Kaufman (1993) for T > 100K, 
and Neufeld, Lepp, & Melnick (1995) for T < 100K. For 
CO and OH, we calculated the cooling rates as in O05 
by using the transition coefficients from Leiden Atomic 
and Molecular Database (Schoier et al. 2005) and fit the 
results by the same form of the formula as in Neufeld & 
Kaufman (1993). The fitting parameters for the CO and 
OH cooling rates are summarized in Appendix A. 



2.2.2. Continuum-emission Cooling 

We solve transfer of continuum radiation by the vari- 
able Eddington factor method with grey approximation 
(Mihalas & Weibel-Mihalas 1984). For a given mean in- 
tensity J, the dust temperature Td is determined by the 
energy balance on the dust grains: 

Kd-B(T d ) = n d J + r co u, (4) 

where the Planck function B(T) = o-bT 4 /it, n d is the 
Planck mean opacity of dust and the heating rate for 
the dust by collisions with gas particles r co n is given by 
(Hollenbach & McKee 1979) 

1 /2 

r coll = 4.4 x 10- 6 (//p) gr n H (— ^) (T - T d ), (5) 

where T is the gas temperature, and (//p) gr is total vol- 
ume of the dust per unit gas mass, which is taken from 
Pollack et al. (1994). For given gas and dust tempera- 
tures, the source function can be written as 

_ _ n d B(T d ) + k s B(T) 



where K g is the gas continuum opacity. Once distribution 
of the source function 5* is given, the mean intensity J 
is found by carrying out standard ray-tracing method. 
We do iterations until the consistent distributions of Td 
and J are obtained. Once radiative transfer is solved, 
the cooling rates by dust and gas continuum emissions 
are given by 

A d = 47r Kd [B(Td)- J] (7) 

and 

A g = 4iTK g [B(T) - J], (8) 

respectively. 

2.2.3. Chemical Cooling and Heating 

The change in chemical binding energy works as heat- 
ing or cooling of the gas. The heating and cooling as- 
sociated the chemical reactions arc treated in the same 
way as in Omukai (2000). We note here that impor- 
tant processes are heating/cooling by the H 2 forma- 
tion/dissociation, respectively. 

2.3. Chemical Processes 

We solve the reduced chemical network of O05, which 
correctly reproduces the evolution of major coolant 
species obtained by a more detailed network calculation. 
Our network includes 24 compounds of hydrogen, deu- 
terium, carbon and oxygen. For the following 15 major 
species, non-equilibrium reactions are solved: H + , e, H, 
H 2 , D+, D, HD, C+, C, CO, C0 2 , O, OH, H 2 0, 2 (see 
Table 1). The other 9 minor species, i.e., H~, CH, CH 2 , 
CO+, 0+, OH+, H 2 0+, H3O+, 0+, which appear only 
as intermediaries of the reactions among the above major 
species, are set to be chemical equilibrium (see Appendix 
B for more details). By performing some experiments us- 
ing the one-zone models, we confirmed that the evolution 
of the major species' abundances in this method agree 
very well with that in the full non-equilibrium treatment. 
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2.4. Initial Core Models 

Initial models for the hydrodynamical calculations are 
prepared using the early evolution in our one-zone model. 
The model is the same as that developed in O05, but 
some rate coefficients of radiative and chemical processes 
are updated to be the same as in the hydrodynamical 
models. Figure iBTl shows the prestellar temperature evo- 
lution calculated by the one-zone model with different 
metallicities. The temperature and chemical abundances 
of the initial models are assumed to be radially constant, 
and their values are taken from the one-zone models at 
10 4 cm~ 3 . For density and velocity distributions, we as- 
sume the critical Bonner-Ebert sphere at rest, which is 
the hydrostatic configuration on the verge of the collapse. 
To initiate gravitational collapse, we enhanced the den- 
sity by twenty percent. The central number density is 
1.2 x 10 4 cm -3 at the beginning of the hydrodynamical 
calculations. In addition to the case of primordial gas, 
Z = 0, we study the cases every one dex in metallic- 
ity from Z = 1CT 6 Z Q ([M/H]=-6) to IZ Q ([M/H]=0). 
The total masses of the dense cores depend on the initial 
temperature at 1.2 x 10 4 cm~ 3 for the hydrodynamical 
calculations; 4.7 ([M/H]=0), 37 (-1), 52 (-2), 1.7 x 10 2 
(-3), 1.6 x 10 3 (-4), 2.3 x 10 3 (-5), and 2.4 x 1O 3 M (-6, 
and Z = 0). The outer boundary is fixed at its initial 
position. 

3. PRESTELLAR COLLAPSE 

Figure IB2I presents the temperature evolution at the 
center of the prestellar cores as a function of the number 
density. The overall evolution is quite similar to that cal- 
culated by the one-zone model (Fig lBll) . justifying the 
one-zone treatment for the core evolution. There are, 
however, small disagreements, in particular, at high den- 
sities and for low metallicity cases. We defer detailed dis- 
cussion on these differences to later in Sec. 3.4, but here 
describe which thermal processes control the tempera- 
ture evolution at each metallicity. The contribution to 
the cooling and heat ing rates by individual processes are 
presented in Figure IB3I for different metallicities. This 
should be compared with Figure 2 of O05, where simi- 
lar plots for the one-zone models are presented. In Fig- 
ure IB51 the effective ratio of specific heat at the center 
7 = dlnp/dhip, which tells the variation of pressure in re- 
sponse to the density variation, is shown for those cases. 
Not e that 7 — 1 equals the gradient of the curve in Figure 
IB2l for constant molecular weight. The effective ratio of 
specific heat is an important index to examine dynamical 
response of self-gravitating clouds to thermal evolution. 
For example, the clouds easily fragments as long as 7 < 1, 
while fragmentation is strongly prohibited for 7 > 1 (Li 
et al. 2003). Another critical value is 7 = 4/3. If 7 
exceeds this value, the dynamical collapse is halted as 
the pressure overcoming the gravity, and a hydrostatic 
object is formed. 

3.1. Thermal Evolution in the Metal-free Case 

In this section, we review thermal evolution of the 
cloud core of a metal-free gas. We then describe the 
effects of metallicity later in Sec. 3.2. We focus on 
deviations from the metal-free case. In the case of 
metallicity [M/H]=-6, metallicity effects are so small 
that the temperature evolution is almost identical to the 



metal-free one except a slight offset at highest densities 
(> 10 20 cm- 3 ). 

Let us summarize here the formation processes of H2 , 
which play a crucial role in the thermal evolution. The 
evolution of H2 concentrations is presented in Figure 
IB41 along with those in the cases with metals. Below 
~ 10 8 cm~ 3 , H2 is formed by the H~ channel: 

H + c^H~+7 (9) 

H~ +H->H 2 +e, (10) 

catalyzed by a small amount of remaining electrons. 
With their recombination proceeding, the H~ channel 
is quenched and the amount of formed H2 saturates at 
~ 10~ 3 (Fig. IB4|) . After this plateau, the H2 abundance 
begins to increase again at ~ 10 8 cm -3 via the three-body 
H2 formation: 

2H + H^H 2 +H, (11) 

and 

2H + H 2 ^H 2 +H 2 . (12) 

All the hydrogen is converted to the molecular form via 
this channel by the density ~ 10 n cm~ 3 . 

Next, let us see the cooling and heating processes (Fig. 
IB3l aV Until very high density ~ 10 19_20 cm~ 3 is reached, 
cooling and heating are always almost balanced, so that 
the evolution is nearly isothermal with temperature dif- 
fering only by a small factor whereas density increases by 
many orders of magnitudes. The effective ratio of spe- 
cific heat 7 remains below 4/3 but is above 1 in this 
period except for brief intervals around 10 9 cm -3 and 
10 11 cm -3 , where 7 falls slightly below unity (Fig. IB5I 
a). The heating is owing to the compression, but for 
10 9 — 10 12 cm~ 3 , where the H2 formation heating associ- 
ated with the three-body reaction (eq. [11] below) domi- 
nates. For the cooling, the H2-line emission contributes 
most until ~ 10 13 cm~ 3 , although some lines become op- 
tically thick at ~ 10 11 cm -3 and this suppresses the cool- 
ing rate gradually towards a higher density. The steep 
decline of the H2 line-cooling rate at 10 16 cm -3 is due to 
the H 2 collision-induced continuum absorption. Another 
molecular species in the metal-free gas, HD, is known 
to play an important role in cooling if a metal-free gas 
is once ionized (Uehara & Inutsuka 2000; Nagakura & 
Omukai 2006; Greif & Bromm 2006; Yoshida, Omukai, 
& Hernquist 2007; McGreer & Bryan 2008). In our case, 
however, it only contributes comparably to H2 at a brief 
period at ~ 10 4 cm~ 3 . 

With gradual increase of temperature, the balance of 
chemical equilibrium between the H2 formation (eq. Ill[) 
and its inverse dissociation reactions shifts in favor of dis- 
sociation in the range 10 13 — 10 15 cm~ 3 . The cooling by 
this subtle dissociation is enough to compensate the com- 
pressional heating, because of the large binding energy 
of molecular hydrogen, \ = 4.48 eV. At ~ 10 15 cm~ 3 , 
H2 formation is somewhat resumed owing to the cooling 
by the H2 collision-induced continuum emission, which is 
a molecular analogue of free-free emission. At densities 
greater than ~ 10 16 cm~ 3 , the cooling rate by this pro- 
cess declines rapidly when the central part of the core be- 
comes optically thick to this same continuum. This time, 
with increase of temperature, the H 2 begins to dissoci- 
ate. The dissociation cooling balances with the compres- 
sional heating until ~ 10 20 cm~ 3 , where the dissociation 
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is almost completed. For higher densities, without effi- 
cient cooling process, the temperature increases almost 
adiabatically and 7 exceeds 4/3. This marks the end of 
the prestellar collapse and the formation of a hydrostatic 
object, i.e. a protostar, at the center (see Sec. 3.3). 

3.2. Metallicity Effects on Thermal Evolution 

Addition of metals alters the thermal evolution of star- 
forming cloud cores significantly in the following ways. 

3.2.1. Atomic and Molecular Metal-line Emission 

For a metal-free cloud core, the efficient radiative cool- 
ing is solely by H2, i.e. its rovibrational line emission and 
collision-induced emission, with some contribution from 
HD around ~ 10 4 cm -3 . With metal enrichment, the 
cooling by line emission of metallic atoms ([CII], [CI], 
[01]) and molecules (CO, OH, H2O) affects the thermal 
evolution. 

Oxygen-bearing molecules, OH and H2O, act as effi- 
cient coolants. In the case of the metallicity [M/H]=-5, 
since the gas temperature remains high at > 200K, OH 
molecules formed at ~ 10 7 cm -3 are readily converted to 
H2O. Cooling by H2O rotational transition exceeds H2 
cooling in the range 10 8 — 10 10 cm~ 3 . With more metals 
of [M/H]=-(3-4), the lower temperature conditions delay 
the conversion of OH to H2O molecules. The line emis- 
sion by OH molecules dominates cooling at 10 6 — 10 8 cm -3 
before the ^-formation heating raises the temperature 
and promotes the conversion to H2O molecules. Af- 
ter that, H2O significantly contributes to cooling until 
10 10 cm- 3 (10 9 cm- 3 ) at [M/H]=-4 ([M/H]=-3, respec- 
tively), when the dust thermal emission becomes impor- 
tant. With higher metallicity, even lower temperature 
due to the CO cooling makes OH and H 2 hard to be 
excited and their importance in cooling diminishes al- 
though, for [M/H]=-2, OH cooling is still important for 
10 6 - 10 7 citT 3 . For [M/H] > -2, the fine-structure line 
transitions of carbon and oxygen are important for cool- 
ing in a diffuse phase (< 10 3 cm~ 3 ), whic h is below our 
initial densities and not shown in Figure IB3I At these 
metallicities, CO becomes the dominant coolant before 
it is overtaken by dust thermal emission. For [M/H]=0, 
this shift happens ~ 10 3 cm~ 3 and the CO cooling does 
not appear in Figure lB3l 

3.2.2. H2 Formation on Dust Grains 

In a metal-free gas, only a small amount of H2 with 
concentration y(H.2) ~ 10 -3 is formed via the gas-phase 
H~ reaction channel before the onset of the three-body 
H2 formation at ~ 10 8 cm~ 3 . In a metal-enriched gas, 
dust grains are usually present as well; efficient H2 for- 
mation proc eeds on the grain surface even at low densi- 
ties. Figure lB4l shows that, at metallicity [M/H] > —5, 
the H2 concentration continues increasing owing to the 
formation reaction on the dust, unlike the metal- free case 
where there is a plateau at ~ 10~ 3 . For higher metal- 
licities, i.e., more dust grains, the H2 formation is more 
efficient, and its concentration is higher at the same den- 
sity before becoming fully molecular by the three-body 
reaction (< 10 10 cm~ 3 ). In the case of [M/H] =-5, for 
example, about three times more H2 are formed by the 
beginning of three-body reaction (10 8_9 cm~ 3 ) than in 
the metal-free case. This higher H2 abundance causes a 



slight deviation in temperature from the metal-free case 
< 10 7 cm -3 even before the small temperature dec reas e 
by the H 2 cooling in the range 10 7 ~ 10 cm~ 3 (Fig. IB2|) . 
With higher metallicity of [M/H] =-(3-4), further H 2 en- 
hancement makes the temperature fall below the thresh- 
old for HD formation at ~ 150K. The temperature de- 
creases by HD line cooli ng un til its rotational levels reach 
LTE at ~ 10 5 cm -3 (Fig. |B2j}. With even higher metallic- 
ity ([M/H] > -2), despite of higher H2 abundance, the line 
cooling by metalic atoms (01) and molecules (CO, OH) 
is more significant than that with H2 and HD molecules. 

Except for the case of the solar metallicity ([M/H]=0), 
H2 formation has not yet been completed by the initial 
density of calculation at 10 4 cm -3 . Heating associated 
with H2 formation has large effects on the temperature 
evolution. In the cases with lowest metallicities ([M/H] < 
-5), the heating by three-body reaction is dominant over 
the compressional heating in the range 10 9 — 10 12 cm~ 3 . 
For metallicity of [M/H]=-4..-l, because of the low tem- 
perature at the onset of rapid H2 formation, the chemi- 
cal heating develops a remarkable bump in the tempera- 
ture evolution, separating two temperature dips, i.e. one 
at lower density by line cooling and the other at higher 
density by dust cooling. In particular, this effect makes 
a large jump in temperature evolution at 10 8 cm~ 3 for 
[M/H] ='-4.. -3. 

It should be noted that the rapid H2 formation phase 
and thus the resultant temperature bump shift towards 
lower density with higher metallicity. Also, the local 
maximum value of the temperature decreases with in- 
creasing metallicity. This is because, for the H2 forma- 
tion in lower density, a larger fraction of the chemical 
binding energy is radiated away before collisions with 
other gas particles make it available for heating, along 
with more efficient radiative cooling by metals. 

3.2.3. Thermal Emission of Dust Grains 

Thermal energy of gas is efficiently transferred via colli- 
sions to dust grains, which has a lower temperature. The 
excess energy deposited on the grain is readily radiated 
as thermal radiation. This works as an efficient cooling 
mechanism at high densities where the collisions bet ween 
gas and dust are frequent enough (see panels in Fig. |B3 ). 
and makes a temperature dip for [M/H] > —5 (Fig. IB2 1. 
The frequent collisions equalize the temperature of gas 
and dust (so-called dust-gas coupling). However, before 
the density increases enough that the dust-gas coupling 
occurs, the dust is cooler than the gas by a large margin. 
Note that we assume there is no external background ra- 
diation. In this case, the cooling rate by dust thermal 
emission obeys Adust = r co n oc Zn^T 3 / 2 (see eq. O, 
where the Z dependence comes from (f /p) g r- The com- 

1 /2 

pressional heating rate is proportional to n H T for the 
collapse at a rate of the free-fall. With increasing density, 
thus, the dust cooling becomes more important relative 
to the compressional heating and finally catches up with 
it. At this point the gas temperature falls down until 
it reaches the dust temperature. Thereafter, the cooling 
rate by dust emission is expressed as Adust oc ZT^ +A in 
the optically thin case, where the Planck mean opacity 
of dust Kd oc . With keeping thermal balance between 
the dust cooling and compressional heating, temperature 

slowly increases as oc ^ ~ n ^ 1 , where (3 ~ 2 is 
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used in the last expression. 

Because of the dependence of Adust on Z, the dust 
cooling becomes less efficient for the lower metallicity 
and then the onset of its efficient cooling is postponed 
until the higher density regime, where the temperature 
is also higher. With metallicity below a threshold value 
([M/H]=-6..-5), the temperature reaches the dust evap- 
oration value before the dust cooling becomes efficient. 
In the case of [M/H]=-6, for example, the dust com- 
pletely sublimates at ~ 1500K (2 x f 14 cm~ 3 ) before its 
cooling reaching the level of compressional heating, al- 
though it does some contribution at 10 13 — 10 14 cm~ 3 . 
The temperature dip by dust cooling is absent in this 
case. With higher metallicities, the dust cooling reaches 
the compressional heating at some point and balances 
it thereafter until the dust thermal emission becomes 
optically thick and ineffective. After the core becomes 
optically thick to the dust absorption, the temperature 
increases almost adiabatically with in creas ing density as 
seen clearly for [M/H] > —4 in Figure lB2l This phase is 
called the first adiabatic phase. We distinguish it from 
another adiabatic phase (second adiabatic phase) at a 
higher density (> 10 21 cm -3 ). Note that a cloud core 
with higher metallicity becomes optically thick to dust 
absorption at progressively lower density and temper- 
ature: for example, 10 n cm~ 3 for Z©, and 10 13 cm -3 
for 10~ 4 Z©. Never theless, evolutionary paths of tem- 
peratures (Fig. IB2j) during and after the first adiabatic 
phase are similar each other for [M/H] > —5 (see Omukai 
2000 for further explanation). Finally, at ~ 1500K, the 
remaining dust species, i.e., silicates and iron, totally 
evaporate. The central part of the core is now free of 
dust and becomes transparent again although the sur- 
rounding matter is still optically thick. Because the sur- 
rounding gas has a lower temperature, the central part 
is still able to cool by continuum emission. This phase, 
however, lasts for a brief interval since the core soon be- 
comes optically thick to the same continuum. The H2 
dissociation begins at this moment and the temperature 
increases only gradually under the action of the dissoci- 
ation cooling. After its completion at ~ 10 21 cm~ 3 , the 
temperature increases adiabatically again. This is the 
second adiabatic phase. The evolution after the begin- 
ning of H2 dissociation is similar to the metal-free case, 
although the temperature is somewhat lower than that 
in the metal-free case. 

3.3. Dynamical Evolution 

We next describe dynamical evolution of the dense 
cores as a response to the thermal evolution seen above 
in Sec. 3.2. In Figure IB61 the rate of enhancement of 
central density p c , which can be given by the inverse of 
the dynamical timescale, t^ n — (dp c /dt)/ p c is shown 
relative to the inverse of the free-fall timescale t^ 1 = 
Vyl^T' By comparing Figures lB5l and lB6l we see an 

anti-correlation between the collapse rate tg/t^ yn and 7: 
where 7 is high (low) , the collapse rate is low (high) ow- 
ing to increasing (decreasing, respectively) pressure with 
collapse relative to the gravity. In particular, when 7 ex- 
ceeds 4/3, the dynamical collapse is halted and the col- 
lapse rate becomes very low with iff/tdyn < 1 (lower dot- 
ted horizontal lines in Fig. IB6|) . indicating that the col- 



lapse is quasi-static, although the core contracts slowly 
by accretion of ambient matter. 

In Figures IB7IIB131 we present radial distributions of 
density, temperature, velocity, and enclosed mass of the 
dense cores at different epochs until the formation of pro- 
tostars for zero-metallicity, and metallicity of [M/H] =-6, 
-5, -4, -3, -2, -1, and 0. Between successive plots, the in- 
crement of central density is 10 25 times that of the pre- 
vious one. The exception is the final phase shown, which 
corresponds to the state slightly after the end of the 
prestellar collapse. Similarity nature of the collapse i s ob- 
viou s in density distributions in Panels (a) of Figurcs [B7I 
IBI31 This collapse closely resembles that of the Larson- 
Penston (LP) similarity solution (Larson 1969; Penston 
1969), which describes the runaway collapse of isother- 
mal clouds. In this solution, a flat-density central part of 
about a Jeans length collapses at a rate of ts/t^ yn = 3.0. 
With increasing density, the Jeans mass and thus the 
matter in the central flat region decreases. Since the col- 

— 1/2 

lapse timescale idyn oc p c becomes shorter and shorter 
with increasing density, a gas in low-density envelope is 
effectively leftover from the evolution of the central re- 
gion during the collapse. The density distribution in the 
envelope has a power of index -2 with 4.4 times that of 
the singular isothermal sphere, while the infall velocity is 
constant at 3.3 times the sound speed in the LP solution. 
The same type of solutions for more general polytropic 
equ atio n of state has been found by Yahil (1980). Fig- 
ure IB6I shows that the collapse rate at the center of our 
cores also becomes close to the value for the LP solu- 
tion, which is iff/idyn = 3.0 for the isothermal case and 
indicated by dotted horizontal line. Note that the value 
iffAdyn higher than unity does not mean the collapse is 
faster than the free-fall, for which tg/td yn — 37r/2 = 4.7 
asymptotically. In some occasions, in particular for the 
cores in the dust-cooling phases, peak values of the col- 
lapse rate exceed the isothermal LP value of 3.0 slightly 
as 7 < 1 in those ranges. For [M/H] =-4 and -3, the tem- 
perature jump by H2 formation heating at ~ 10 8 cm" 3 
is observed as a sharp peak in 7 in Figure lB5l The slow- 
down of the collapse at this epoch re sults in a tempo- 
rary small value of tg/td yn ~ 0.5 (Fig. IB6I as downward 
spikes) , which indicates that the collapse becomes quasi- 
static. Although the sudden slow-down of the colla pse 
dev elops a shock around the heated core (see Figs. IB9I 
and lBlOb ). this heating phase lasts only for a brief period 
and then dynamical collapse is soon recovered. While H2 
formation is in progress, small oscillations in the central 
density are imposed by the associated heating on the 
overall collapse. These oscillations are seen for [M/ H]=- 
1 as notches around the peak of the bump in Figure IB21 
For other cases, they are no t no ticeable in Figure IB21 
but still observable in Figure IB3I as notches in the com- 
pressional heating rates at H2 formation epoch. After 
the core becomes optically thick, temperature increases 
adiabatically and 7 exceeds 4/3. The pressure overcomes 
the gravitational pull before long and a hydrostatic core 
forms. This hydrostatic core is often referred to the first 
protostellar core in literature as it appears before the gen- 
uine protostar, which is also called the second protostel- 
lar core. Although dynamical collapse is halted in this 
phase, the core continues contracting on average, but its 
oscillates acoustically many times, in response to the in- 
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creasing mass of the hydrostatic core by accretion of the 
envelope material. Finally, at ~ 10 16 cm~ 3 , dynamical 
collapse is resumed owing to the cooling by H2 dissoci- 
ation. When the dissociation has almost completed, the 
temperature again begins to increase adiabatically. A 
hydrostatic core formed in this second adiabatic-phase 
grows eventually to a star and thus is called the proto- 
star. The formation of the protostar marks the transition 
from the prestellar phase to the protostellar phase. Ob- 
servationally, the protostars still in the early protostar 
phase are considered be the Class objects, while the 
Class I objects are later but still in the protostar phase. 

In the metal-free case, the temperature increases grad- 
ually at almost a constant rate from 10 4 cm -3 to 
10 19 cm -3 without remarkable features. Consequently, 
the radial density profile does not deviate noticeably from 
that of the LP-type similarity solution until > 10 20 cm" 3 , 
where the appearance of a hydrostatic core results in 
steepening of density gradient. For a polytropic gas 
whose equation of state is P oc p 7 , the envelope den- 
sity profile obeys p oc r~ 2 /( 2 " 7 ) (Larson 1969; Yahil 
1983). Since 7 takes the values of 1 - 1.2 at < 10 19 cm" 3 
(Fig. IB5|) . the density gradient in the envelope is ~ —2.2, 
which is slightly steeper than the value for the isother- 
mal collapse, —2. The three-body H 2 formation makes a 
sma ll dip in temperature evolution at ~ 10 10 cm~ 3 (Fig. 
IB2|) . which is observed also in radial temperature dis- 
tribution at 10 2_3 AU (Fig. lB7l b). The formation of 
the protostar at 10 21_22 cm -3 results in sudden halting 
of the inflow, which develops th e ac cretion shock at the 
stellar surface ~ 10 _2 AU (Fig. IB7I c). Also, the proto- 
star formation can be seen as steepening of density profile 
around the protostellar surface, rapid increase of temper- 
ature, and a plateau in the enclosed mass profile. The 
mass of the protostar at its birth is somewhat smaller 
than ~ 10 -2 Mq. With metal enrichment as small as 
[M/H]=-5, the density, velocity, and mass profiles look 
similar to those in the metal-free case, although tempera- 
ture profile deviates from the metal-free one at ~ 10 3 AU 
due to H 2 cooling and ~ 1AU due to dust cooling (Fig. 
IB8[) . Owing to the difference in temperature evolution 
for > 10 16 cm~ 3 (Fig. IB2|) . the effective adiabatic index 7 
exceeds the critical value 4/3 at 10 21 cm~ 3 for [M/H]=-5 
and higher metallicity cases, while at 10 20 cm -3 for Z = 
(see Fig. IB5|) . Thus the formation of protostar is de- 
layed to slightly higher density and, the mass and radius 
at protostar formation are somewhat smaller than the 
metal-free case. 

For the case of [M/H]=-4, the heating by H2 formation 
make s a j ump in temperature distribution at ~ 10 3 AU 
(Fig. IB9b ). This causes slow-down of collapse and for- 
mation of a short-lived hydrostatic core, whose appear- 
ance can be seen as kinks in the density and mass profiles 
(Fig s. IB9a . d), and a shock in the velocity profile (Fig. 
IB9b ). A temperature dip by dust cooling is obvious at 
~ 10AU. Almost simultaneously with the dust-gas cou- 
pling, the core becomes optically thick to the dust opac- 
ity. Subsequently, the temperature increases rapidl y fo r 
10 14 - 10 16 cm~ 3 , but 7 does not exceed 4/3 (Fig. [B5|) . 
The dynamical collapse continues although slowing down 
with iff/idyn reaching unity at 10 17 cm~ 3 , w hich indicates 
the collapse is now quasi-static (Fig. IB6|) . This slow- 
down is seen as a small spike at ~ 0.1 AU in the velocity 



profile (Fig. IB9b ). Even in a non-spherical collapse, a 
cloud becomes round by this heating and fragmentation 
in the later dust-cooling phase might be suppressed. If 
so, low-mass fragments and stars cannot be formed in the 

range [M/H] 4. Tsuribe & Omukai (2008) attributed 

the origin of apparent paucity of stars in the range [M/H] 
= —(4.8. .4.0) in the halo to it. The evolution after the 
H2 dissociation > 10 16 cm~ 3 is identical to the case of 
[M/H] =-5. 

The case of [M/H] =-3 is qualitatively similar to 
[M/H] =-4 (Fig. iBlO]) . However, the H 2 formation heat- 
ing, as well as the temperature increase in the opti- 
cally thick regime, have larger effects. Even before the 
rapid H 2 -formation heating at 10 8 cm~ 3 , slower heating 
has started and temperature increases gradually with 
7 = 1.1 — 1.2, higher than ~ 1 for [M/H]=-4 in the same 
regime. The collapse then is slower with tff/idyu — 1-5 
while it is ~ 2 for [M/H] =-4. Thus the consequence 
of the rapid heating is almost complete halt of the col- 
lapse for a brief period, as observed in a deep trough in 
the velocity profile at 2000AU (Fig. lBl~0b ). The central 
part becomes optically thick to the dust absorption at 
10 12 cm~ 3 , and the adiabatic contraction phase continues 
until 10 16-17 cm~ 3 , where the dust sublimates and then 
the H2 dissociation follows. Note that since the adiabatic 
phase begins at lower density for higher metallicity and 
ceases at almost the same density, the ad iabatic phase 
is more prolonged. As seen in Figure IB61 the collapse 
becomes quasi-static for > 10 15 cm -3 and a hydrostatic 
object - the first protostellar core - forms. Signs of the 
hydrostatic core formation, i.e., density steepening, tem- 
perature increase, velocity dip, and flatteni ng of enclosed 
mass profile, are seen at ~ 1 AU i n Fig ure IB10I 

In the case of [M/H] =-2 (Fig. [Bill , the H 2 forma- 
tion heating is gradual with 7 ~ 1.3 without a jump in 
temperature. Thus, no hydrostatic object forms in this 
phase. The hydrostatic objects (first hydrostatic core) 
forms at 10 14 cm~ 4 , after the central part becomes opti- 
cally thick to dust absorption. The core starts collapsing 
dynamically again by dust evaporation and H 2 dissocia- 
tion. Finally, at 10 22 cm -3 , the second hydrostatic core, 
i.e. protostar, forms with the same condition as in other 
cases with [M/H]> -5. 

For higher metallicities (Figs. IB12l and lB13|) . the qual- 
itative nature of the evolution remains the same, but 
the onsets of both the efficient dust-cooling and adia- 
batic phases are shifted to lower densities. Because of 
the longer adiabatic phase, the appearance of the first 
hydrostatic cores become more conspicuous (at ~ 1AU). 
The radius of the first protostellar core becomes slightly 
larger with increasing metallicity. Note that for the so- 
lar metallicity ([M/H]=0) case, the initial (at 10 4 cm~ 3 ) 
temperature is already low, t he ou ter bump is missing in 
the temperature profile (Fig. IB13b ). 

3.4. Comparison with the One-zone Model 

Figures IB1I and IB2I show some differences between the 
hydrodynamical and one-zone models, despite the overall 
similarity. In general, temperatures tend to be higher in 
the hydrodynamical models than in the one-zone models. 
The differences are summarized as follows; (1) In the one- 
zone model, all the evolutionary tracks converge for > 
10 16 cm -3 , while in the hydrodynamical models, the two 
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lowest metallicity cases (Z = and [M/H]=-6) remain 
higher than in others. (2) At [M/H]=-5, the temperature 
dip by dust cooling around 10 14 cm -3 in hydrodynamical 
model is not as obvious as in the one-zone model. (3) 
At [M/H]=-3, H2 formation heating results in a large 
temperature jump at ~ 10 8 cm~ 3 in the hydrodynamical 
model, while in the one-zone model it ends up with a 
gradual increase in temperature. 

All these deviations can be attributed to the fact that 
the collapse rate of the core tends to be faster in the 
hydrodynamical model than that assumed in the one- 
zone model. In the one-zone model, the dynamical 
timescale is assumed to be equal to the free-fall time, 
i.e., tg/tdyn — 1. In the hydrodynamical models, soon 
after the beginning of collapse from the rest at 10 4 cm~ 3 , 
the collapse rate exceeds ts/td yn > 1, which indicates 
that the collapse is highly dynamical. It should be noted 
that the free-fall timescale is the time needed for the den- 
sity of the homogeneous gas at rest to reach infinite if the 
pressure is negligible, and iff/idyn higher than unity does 
not mean that the collapse is indeed faster than the free 
fall. 

Observing convergence of temperature tracks in high 
densities in the one-zone model, Omukai (2000) demon- 
strated analytically that the cores that cool by contin- 
uum emission and collapse over a timescale proportional 
to the free-fall time have the same specific entropy when 
they become optically thick to the same continuum, re- 
gardless of their different metallicities. The same value 
of specific entropy assures the convergence of the evo- 
lutionary tracks after the cores become optically thick. 
In the discussion of Omukai (2000), however, the gas was 
assumed to be fully molecular when the core becomes op- 
tically thick to the continuum. In our hydrodynamical 
models, due to higher temperature than in the one-zone 
models for Z=0 and [M/H]=-6, about 10% of hydrogen 
molecules are already dissociated when the core becomes 
optically thick to the H2 collision-induced absorption. 
The lower binding energy and thus higher specific en- 
tropy result in higher temperatures even after the core 
becomes optically thick. 

4. EVOLUTION OF PROTOSTARS IN THE ACCRETION 
PHASE 

In Sec. 3, we have seen the collapse of dense cores 
and formation of tiny protostars at the center. Although 
very small mass of < 10 _2 M Q at its birth, the protostar 
is surrounded by a large amount of gas in the envelope, 
which is a remnant of the parental core. The protostar 
thus grows in mass by orders of magnitude through accre- 
tion of these materials. This is the main accretion phase, 
corresponding to Class and I phases, observationally. 

When the protostar is not very massive (< 10 — 
IOOMq), stellar feedback to the accretion flow is not sig- 
nificant yet. Then the accretion rate is determined by 
the density and velocity structures in the envelope at 
the moment of protostar formation, which are shown in 
Figures TB14I and IB151 respectively. For density distribu- 
tions, the number density multiplied by r 2 is presented 
to emphasize differences among the cases with different 
metallicities. As a general trend, both the density and 
velocity are higher for lower metallicity owing to higher 
prestellar temperature as is expected for the LP simi- 
larity solution, which has density 4.4 times the singular- 



isothermal-sphere value cx c 2 /r 2 tx T/r 2 and infall veloc- 
ity 3.3 times the sound speed c s oc T 1 / 2 (Sec. 3.3). We 
evaluate the accretion rate as 



M{M r ) = 4irr 2 (-v)p 



(13) 



M,. 



These values are shown in Figure IB16I as a function of 
the enclosed mass M r , which can be regarded as the pro- 
tostellar mass M* at the time when the mass element 
of interest falls on the protostar. It s hould be noted 
that the accretion rates shown in Figure IB16I are about 
ten times higher than the often cited values of c 3 /G, 
which are valid for the accretion from the static singular 
isothermal sphere (Shu 1977), because of the dynamical 
nature of the collapse. For the LP solution, this factor 
is ~ 50. Although the density distribution soon reaches 
that of LP solution, convergence of the velocity to the LP 
solution is slow and does not reach the value of LP solu- 
tion during the actual collapse. This results in a smaller 
factor of ~ 10 enhancement instead of ~ 50 from c 3 /G 
(Foster & Chevalier 1983; Larson 2003) 

In the case of metallicity [M/H]=-3, expansion, rather 
than infall, is observed at a small radial range owing to 
oscillation of the core at the sudden H2-formation heat- 
ing. This results in the negative value of accretion rate 
by equation (fl3l thin line in Figure [BTBl at ~ 10M Q ) and 
followed by enhanced rate later. In this range, the rate 
does not give the correct accretion rate. In the follow- 
ing calculation of protostellar evolution, we use a rate 
extrapolated from lower protostellar mass until the ac- 
cretion r ate be comes the same as this value (dotted line 
111 Figure [B16l) in this range. 

With the accretion rate obtained above, we calculate 
the accretion evolution of protostars as in Hosokawa & 
Omukai (2009a). In this method, originally developed 
by Stahler et al. (1980) and Stahler & Palla (1986), a 
protostar model is constructed by solving ordinary stel- 
lar structure equations assuming the hydrostatic equi- 
librium, while the envelope is solved under the assump- 
tion of the stationary accretion with a given accretion 
rate (s ee Ho sokawa & Omukai 2009a for more details). 
Figure IB 1 71 shows the evolution of protostellar radius 
with increasing protostellar mass for different metallic- 
ities. The calculations are initiated with the protostar 
mass of 0.01M Q . For [M/H]<-5, the accretion rate is 
too high at the beginning and we cannot find conver- 
gent solutions for protostars by our numerical scheme. 
Then we start from a smaller rate and increases it wit h 
protostellar mass until the rate shown in Figure IB16I is 
reached. For those cases, we show the protostellar radii 
only after the correct accretion rates are reached and 
transient behaviours disappeared. Since the accretion 
rates adopted are similar, this result resembles closely 
that in Figure 3 of Hosokawa & Omukai (2009b), which 
shows the cases of 10c 3 /G, calculated by t emperatures 
from the one-zone model. From Figure [B 171 we first no- 
tice that lower metallicity protostars have larger radii 
owing to higher accretion rates. Since a short accretion 
timescale does not allow accreted matter to cool enough 
near the stellar surface before embedded inside, a proto- 
star with higher accretion rate has higher entropy con- 
tent and thus larger radius. In all cases in Figure IB17[ 
the protostellar radius initially increases gradually with 
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addition of newly accreted matter. When the protostar 
is small in mass, the opacity due to the free-free absorp- 
tion is so high inside the star that the embedded mat- 
ter has no time to release its heat (so-called the adia- 
batic accretion phase). The protostar with solar metal- 
licity ([M/H]=0) has small entropy content and thus the 
entropy generation by deuterium burning influences the 
evolution, which is observed as an enhanced increase rate 
of radius at ~ 0.4M Q . In other cases, the deuterium 
burning has no large impact in their evolution. The grad- 
ual increase in radius in the adiabatic accretion phase is 
followed by abrupt swelling. With increase of protostel- 
lar mass, the temperature inside also increases and the 
opacity decreases as oc T~ 35 . At some point, the de- 
creased opacity allows the heat kept in the stellar interior 
to propagate outward by way of radiative diffusion. At 
this moment, the protostellar radius swells abruptly and 
thereafter decreases with throwing away the contained 
heat (so-called the Kelvin- Helmholtz (KH) contraction). 
Finally, when the heat generated by hydrogen burning 
compensates the radiative loss, the star re aches stable 
hydrogen burning phase (solid circles in Fig. IB17|) of the 
zero-age main sequence. If the accretion is not halted by 
other mechanism, for example, blow away of the envelope 
by radiative force onto the dusty envelope, quenching of 
mass reservoir in the envelope, etc, the protostar con- 
tinues to grow in mass. The radius after the arrival of 
zero-age main sequence is larger for the higher metallic- 
ity star since more efficient hydrogen burning by the CN 
cycle halts the KH contraction earlier. 

For cases with [M/H]< —5, the total luminosity, i.e. 
the sum of the contributions from the accretion L acc 
and interior one Lj nt , becomes close to the Eddington 
luminosity during the KH contraction because of the 
high accretion rate and decreasing radius; note that 
iacc/^Edd oc M/R*. At this moment, the surface layer 
of the protostar is thrust outward and the protostellar ra- 
dius begins to inflate violently. This marks the end of the 
stationary accretion. The accretion may stop completely 
if this eruption clears away all the matter in the enve- 
lope. Another possibility is that the mass accretion is 
regulated by this effect. If the accretion rate is temporar- 
ily decreased, the luminosity falls below the Eddington 
value. Thus the accretion may continue in a sporadic 
fashion. For other cases, we continue calculation until 
all the matter in the parental core is accreted. The ar- 
rival to the main-sequence occurs at higher protostellar 
mass for lower metallicity, e. g., at 5OM for [M/H]=- 
4, while at 4M© for [M/H]=0. For [M/H]=0, due to a 
small amount of the mass reservoir (4.7 Mq), the accre- 
tion ceases shortly after the arrival to the main-sequence. 

5. DISCUSSION 

5.1. Radiative Transfer Effects 

In this paper, we have solved radiative transfer for the 
continuum by ray tracing method, while in the one-zone 
model we just reduced the cooling rate from the optically 
thin value A COIrt)t hm by 

Acont = A contithin min(l,T~ 2 ), (14) 

where r is the optical depth of the continuum absorp- 
tion. To see whether such a simple and local treatment 
works in hydrodynamical models, we also calculate the 



cases where the continuum cooling rate at radius r is 
evaluated by using equation (fT4|) . where outward radial 
optical depth r r is used for r. The evolution of cen- 
tral temperature is shown in Figure IB 181 (dashed lines), 
along with those with radiative transfer (solid lines) for 
comparison. The central evolution is reproduced very 
well by this simple treatment although the temperature 
is slightly lower near the opacity limit r ~ 1 and higher 
when dust is evaporated at ~ 10 16 ~ 17 cm~ 3 . The former 
is likely due to our artificial switching of the cooling rate 
by equation (fT4|l from optically thin to thick cases. In 
the radiative transfer cases, this transition proceeds in 
a more gradual way since the cooling rate is somewhat 
reduced from the thin value for the optical depth < 1. 
For the Z — case, higher cooling rate near the opacity 
limit at ~ 10 16 cm -3 results in a slightly higher molec- 
ular fraction, which mitigates temperature increase in a 
later phase and delays the protostar formation (see Fig. 
IB19[) . The higher temperature around the dust evapo- 
ration is explained as follows. In the dust evaporation 
phase, the central part is already devoid of dust and 
transparent, while the surrounding matter is still opti- 
cally thick. Then the central part can deposit its heat to 
the cooler surroundings by the gas continuum radiation. 
In this case, the rate by equation (TT4"]) under-estimates 
the real value since the radial optical depth r r , which in- 
tegrates all the opacity along the radial direction, is too 
large. This results in the higher temperature. 

Our ray-tracing method works well to describe the 
core evolution accurately. However, it does not to re- 
produce radial temperature distribution very well. The 
radial temperature distribution at the protostar forma- 
tion is presented in Figure [Bl~9l (dashed lines) for Z = 0, 
[M/H]=-4, -2, and 0. We see clearly deviations in tem- 
perature near the first protostellar core at 1-100AU. In 
the radiative transfer model, the temperature are higher 
owing to radiative heating by the inner warm gas. Higher 
entropy of this heated gas results in larger size of the first 
protostellar core in a later phase than evaluated by the 
central temperature (e.g., Tomida et al. 2010). 

5.2. Critical Metallicity for Dust-induced Fragmentation 

Armed with the temperature evolution obtained above, 
we now discuss the critical metallicity for fragmentation 
to produce low-mass objects. Tsuribe & Omukai (2006) 
studied the fragmentation of dense cloud cores in the 
dust-cooling phase by three-dimensional numerical sim- 
ulations. They followed the evolution of the core elonga- 
tion £ = a/b — 1, where a (b) is the semi-major (minor) 
axis. In their simulation, a core with £ ~ 1 at the onset 
of the dust-cooling phase fragmented to multiple clumps 
when £ reached a critical value £ ~ 30, which corre- 
sponds to ~ 3 in the linear theory. For more details, we 
refer Tsuribe & Omukai (2006). We apply this criterion 
to our hydrodynamic simulation as follows. First, we as- 
sume that the cores have modest elongations of about 
unity at the beginning of the dust-cooling phase (defined 
where 7 < 1 and cooling is dominated by the dust). We 
expect that the core fragments in this phase when the 
elonga tion e xtrapolated by the linear theory exceeds 3. 
Figure [B20l shows the elongation factor, which is defined 
as how many times the elongation has grown from its 
initial value for [M/H]=-l...-5. For [M/H]=-6, there is 
no clear dust-cooling phase, while for [M/H]=0, rapid 
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dust-cooling phase with 7 < 1 has almost finished by the 
beginning of the hydrodynamical calculation at 10 4 cm~ 3 . 
For cases [M/H]> —5, the elongation factor well exceeds 
the critical value of 3, while for the case with [M/H]=-5, 
it only marginally reaches 3. For lower metallicities, it is 
clear that the elongation factor remains below 3. We thus 
conclude that for cores initially elongated with £ ~ 1, 
the critical metallicity is [M/H] =-5 in our hydrodynam- 
ical models. Of course, whether fragmentation actually 
occurs depends on our assumption that the elongation 
is of order of unity at the beginning of the dust-cooling 
phase. In addition, in the discussion above, we used the 
temperature evolution derived by the calculation assum- 
ing spherical symmetry. When the core becomes signif- 
icantly elongated, the collapse timescale as well as the 
cooling timescale will be different, which alters the tem- 
perature evolution. Note also that the dust-cooling rate 
depends on the uncertain dust-gas ratio and dust compo- 
sition in the early universe (Schneider et al. 2006). For 
example, if the dust-gas ratio is lower than assumed here, 
the critical metallicity should be shifted upward with its 
inverse proportion. Also more refractory dust species, 
e.g., amorphous carbon rather than organic carbon, can 
survive in higher temperature condition of lower metal- 
licity gas, thus enhancing the dust-cooling rate. Despite 
these uncertainties, the value of critical metallicity de- 
rived here provides an approximate estimate for further 
studies of low-metallicity star formation. 

5.3. Protostellar Heating to Prevent Fragmentation 

Stellar radiation heats up the surrounding material. 
For example, Vazquez-Semadeni et al. (2010) studied 
numerically the effect of the stellar ionization feedback 
from massive stars on the evolution of giant molecular 
clouds, and found that it controls their star formation 
efficiency. Since the Jeans mass gives the typical mass of 
fragments, this effect suppresses production of many tiny 
objects as a result of fragmentation. In the context of for- 
mation of present-day massive stars, Krumholtz (2006) 
pointed out that heating by accretion luminosity from a 
protostar alleviates fragmentation of the parental core, 
which facilitates monolithic collapse to a single massive 
object rather than fragmentation into a cluster of small 
ones. Krumholtz et al. (2007) confirmed this effect by ra- 
ditation hydrodynamical calculation. Also, Bate (2010) 
found that radiative feedback from an accreting proto- 
star suppresses fragmentation of the circumstellar disk 
by radiation hydrodynamical simulation. Such a radi- 
tive feedback could be potentially important also in a 
low-metallicity star forming core, owing to high accre- 
tion rate and thus luminosity of the central protostar 
and close separations between the fragments in the dust- 
cooling phase. Here we examine whether this effect sup- 
presses fragmentation in low-metallicity environments, as 
studied by Krumholz (2006) for the present-day star for- 
mation. 

In the envelope of a protostar radiating at luminos- 
ity L, the dust is heated to a temperature Th ea t = 

(L/167rcrr 2 ) , which comes from the condition of ra- 
diative equilibrium. Note that this is just a measure 
for heating effect and is neither the actual dust nor gas 
temperatures: the true values of those temperatures are 
determined by taking account of the collisional coupling 



between the gas and dust. In the following, we com- 
pare the gas temperature T env in the envelope before the 
heating and Thcat, which is the maximum temperature 
by the heating effect in the region where Th ea t > T env . 
This dust temperature Thcat is shown in Figure IB21I for 
log(M*/M ) = -1.5,-0.5,0.5 and 1.5. The innermost 
radii of the lines for Th ca t correspond to the positions of 
rarefaction waves. During the accretion phase, the inner 
matter falls onto the protostar earlier. The rarefaction 
front propagates outward with a sound speed and locates 
at ~ 2M„ in terms of the Lagrangian mass coordinate at 
the time of protostellar mass for the density distribu- 
tion p cx r~ 2 . The matter inside of the rarefaction front 
is falling at about a free-fall speed and probably is not 
affected by the protostellar heating anymore. In Figure 
IB2U we see that Theat at a fixed radius increases with 
the protostellar mass due to increasing luminosity. How- 
ever, the propagation of the rarefaction front is faster 
than increase of heating effect: the temperature at the 
rarefaction front decreases with increasing protostellar 
mass. Also, the effect of heating is less pronounced for 
lower-metallicity protostars. Below [M/H] < —3, Th ca t 
falls below T env at anytime outside the rarefaction front. 
Namely, the protostellar heating has no impact on the 
evolution because of the already high prestellar temper- 
ature. For [M/H] =-2 and -1, the radiative heating may 
fill the temperature dip due to the dust cooling, but this 
phase lasts only for a short period. By = O.3M , this 
radiative heating effect becomes negligibly small. Only 
for [M/H]=0, the temperature due to the radiative heat- 
ing is higher than the prestellar value throughout the 
evolution. 

The tendency that radiative heating has more impact 
for higher metallicity (i.e., lower prestellar temperature) 
and for lower protostellar mass can be understood in the 
following way. When a protostar is not massive and still 
in the adiabatic phase, i.e., < 1OM , its radius is 
approximated as (Stahler & Palla 1986) 



R* = 48. IR 



© 



0.27 



Al 



4.41 x 10- 3 M Q /yr 



We write the accretion rate as 

m=4 



(15) 



(16) 



= 1.63 x 10- 5 A/ Q /yr 1O T^ 2 , (17) 

by using a numerical factor <j) (~ 10 in our case), and 
0io = 0/10 and T10 = T/10K. In the adiabatic phase, 
the protostellar luminosity is dominated by the contri- 
bution from the accretion luminosity: 



GM*M 
i?* 



(18) 



= l.OlxlO 2 T 0? o 59 T 1 o o 89 ™«- 73 , (19) 

where to* = M*/M . Thus, temperature of dust heated 
by the protostar is 



1/4 



= 2.82 x 10 2 K 0V 5 T° 22 m°- 18 rr 1/2 , (21) 
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where r±o — r/lOAU. Suppose the location of rarefaction 
front M rw = 2M* and density distribution is a factor of 
a times that of the singular isothermal sphere p = apsis ■ 
For the Larson-Penston collapse a ~ 4.4, and we set it a 
representative value: 0:4.4 = a/4.4. Using the radius of 
the rarefaction front 



CtCg 

= 5.59 x 10 3 AU all T io m *- ( 23 ) 

We obtain the dust temperature due to protostellar heat- 
ing at the rarefaction front 

T hcat (r rw ) = 11.9K all ^ T? 72 m^ 32 . (24) 

Therefore, the condition on the protostellar mass for 
?heat(7rw) to exceed the prestellar temperature is given 

by 

M. < 1.7M alf ^ 46 T 10 - 88 . (25) 

This demonstrates that the radiative heating effect is im- 
portant only up to some mass given by equation (|25j) . 
which decreases rapidly with increasing prestellar tem- 
perature, as has been discussed from Figure IB211 

In the course of the accretion, circumstellar disks form 
owing to the angular momentum. The gravitational in- 
stability of rapidly accreting disks were investigated by 
Kratter et al. (2010) for the isothermal gas by way of 
global numerical experiments. They found that the cen- 
tral parameter for the stability of an accreting disk is 
£ = M /(c 3 a/G) = (/>(c s /c Si d) 3 , which parameterized the 
accretion rate. Here c Si d is the sound speed in the disk 
and c s is that in the parental core, i.e., before the matter 
falls on the disk. Rapidly accreting disks with £ > 3 be- 
come unstable and fragment into multiples or binaries, 
although the exact value for threshold depends on a ro- 
tational parameter of the parental cores, which stabilizes 
the disks somewhat. In our cases, since (f> ~ 10, the disks 
will become unstable and fragment without significant 
heating. Being closer to the star, cicumstellar disks can 
be more vulnerable to the radiative feedback than the 
matter outside the rarefaction wave, which has been dis- 
cussed above. Although not discussed here, the radiative 
feedback onto the disk is interesting in this respect and 
awaits further investigations in future. 

6. SUMMARY 

We have studied star formation with different metallic- 
ities, — IZq, under the assumption of spherical symme- 
try. Until the birth of hydrostatic protostars, the gravi- 
tational collapse of a dense cloud core is calculated using 
chemo-radiation hydrodynamics. The subsequent proto- 
stellar accretion is followed by assuming the accretion 
flow is stationary while ordinary stellar structure equa- 
tions are solved for the central protostar. 

Our findings are summarized as follows: 

(1) Prestellar temperature is higher in a lower- 
metallicity cloud core during optically thin stage, while 
it approaches a single density-temperature relation after 
the core becomes optically thick. Only in lowest metallic- 
ity (< 10~ 6 Zq) cases, the temperature remains slightly 
higher than other cases. 

(2) Since the temperature evolution converges in cases 
with metallicity > 1O~ 5 Z , physical conditions of pro- 



tostars at their formation are universal with a mass of 
10~ 3 Af Q , a radius of O.5i? . The density and tem- 
peratures are ~ 10 22 cm~ 3 and 3 x 10 4 K respectively. 
In the cases with < 10~ 6 Z@, owing to the slightly 
higher prestellar temperature, the protostar formation 
occurs about an order of magnitude earlier in density at 
~ 10 21 cm -3 , with slightly lower temperature 2 x 10 4 K, 
and somewhat higher mass 2 x 1O _3 M and radius IRq. 

(3) The prestellar collapse approaches the Larson- 
Penston self-similar solution, for which the velocity is 
3.3 times the sound speed c s and the density is 4.4 times 
the singular isothermal sphere value (oc c 2 ). Reflecting 
the higher prestellar temperature, the velocity and the 
density in the envelope of a lower-metallicity dense core 
is higher at a given radius. 

(4) The temperature evolution at the cloud center is 
roughly consistent with that found in our previous one- 
zone models, although collapse is faster than assumed 
in the one-zone model, which results in a slightly higher 
temperature at the center. The faster collapse, in partic- 
ular, affects the temperature in the dust-cooling phase. 
Since the fragmentation into low-mass objects occurs in 
this phase, the critical value of metallicity for inducing 
the fragmentation can be altered by this effect. 

(5) The accretion rate after the birth of the protostar 
is higher in the lower-metallicity case because of the re- 
lation between accretion rate M and the prestellar sound 
speed c s in similarity solutions, M = c/>c 3 /G, where <f> is 
a non-dimensional constant. For the isothermal Larson- 
Penston solution, which describes a highly dynamical col- 
lapse, (j> — 46.9, while for the Shu solution, which is valid 
for accretion from a static isothermal sphere, <f> = 0.975. 
In our cases, owing to modestly dynamical nature of the 
collapse, the accretion rates fall in between these two 
extremes and ~ 10c 3 /G. 

(6) The evolution of protostars under the accretion 
rates obtained from the prestellar collapse calculations 
is followed under the assumption of stationary accretion 
flows. The protostellar mass at its main sequence ar- 
rival is higher in the lower-metallicity case owing to the 
higher accretion rate: for example, it is 4M for 1Z©, 
while 50M© for 10~ 4 Z Q . For < 1O~ 5 Z , too high an ac- 
cretion rate causes the luminosity to reach the Eddington 
limit and the accretion is almost halted at 80M Q . No sta- 
tionary solution exists thereafter. The accretion may be 
halted temporarily. Once the accretion rate is reduced, 
however, the luminosity falls below the Eddington limit 
and thus the accretion would be resumed. The accre- 
tion may continue in such a non-stationary or a sporadic 
fashion. 

(7) The critical metallicity for low-mass fragmenta- 
tion is evaluated by applying the criterion by Tsuribe & 
Omukai (2004) to the calculated temperature evolution 
at the center. With metallicity Z > 1O~ 5 Z , a mod- 
erately elongated cloud core is likely to fragment into 
low-mass objects during the dust-cooling phase. 

(8) Radiation from a forming protostar heats up its 
ambient medium, thereby possibly preventing the frag- 
mentation in the surrounding region. This protostellar 
heating effect is evaluated for different metallicity cases. 
For low-metallicities < 1O _2 Z , the prestellar tempera- 
ture is already high, and thus the heating effect is turned 
out to be insignificant. For higher metallicities, proto- 
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stellar heating can be important especially in the early 
protostellar phase. 

Effects such as rotation, magnetic fields, and turbu- 
lence, which are prevalent in the nearby star-forming 
clouds, have not been addressed in our spherical sym- 
metric analysis. They would delay the contraction dur- 
ing the collapse phase and modify the accretion rate in 
the subsequent phase. In addition, rotation facilitates 
the fragmentation by forming a thin disk, while mag- 
netic fields counteract it by extracting the angular mo- 
mentum (Machida et al. 2008; Hcnnebelle & Teyssier 
2008). Those effects surely await future studies. 

Overall, our calculations clarified the similarlity and 
differences in the formation of protostars with different 
metallicities. While some of the important issues, such 
as the critical metallicity for fragmentation, need further 



studies using full three-dimensional simulations, our re- 
sults provide a foundation for such subsequent studies. 
Low-metallicity star formation is particularly important 
in the context of the formation of premeval galaxies. 
Most of the planned ground-based and space-borne tele- 
scopes are aimed at observing the first galaxies in which 
the gas metallicity is supposed to be very low. Theo- 
retical studies such as those presented here will give an 
invaluable input for the future observations. 
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APPENDIX 

A. REDUCED CHEMICAL NETWORK 

The set of reactions are presented in Table 1. The number of reactions and rate coefficients are the same as Table 
1 of Omukai et al. (2005). This model includes 24 compounds of hydrogen, deuterium, carbon and oxygen, i.e., H+, 
e, H, H2, D + , D, HD, C + , C, CO, CO2, O, OH, H9O, O2 H~. Non-equilibrium reactions are solved for the following 
15 major species: H+, e, H, H 2 , D+, D, HD, C+, C, CO, C0 2 , O, OH, H 2 0, 2 (see Table 1). The other 9 species 
(H~, CH, CH 2 , CO + , + , OH+, H 2 + , H 3 + , Oj ), which appears only as intermediaries of the reactions among the 



above major species, are set to be chemical equilibrium as follows: 

H : k 2 n(H)n(e) = k 3 n(H~)n(H) (Al) 

CH : [k 37 n(H) + k 2H n(H 2 )]n{C) + k 22 n(CH 2 )n(H) = n(C H)[k 21 n{H) + k 31 n(H 2 ) + k 42 n(0)} (A2) 

CH 2 : k 62 n(H 2 )n(C) + k 31 n(H 2 )n{C H) = k 22 n(CH 2 )n(H) (A3) 

CO+ : k 48 n(C + )n(OH) = k 54 n(CO + )n(H) (A4) 

0+ : k 30 n{H + )n(O) + k 49 n(C+)n(0 2 ) = [k 50 n{H) + k 51 n(H 2 )}n(0+) (A5) 

OH+ : k 51 n(0 + )n{H 2 ) + k 45 n{H+)n(OH) = [k 52 n(H 2 ) + k m n(e)]n{OH+) (A6) 

H 2 0+ : k 52 n(OH+)n(H 2 ) + k 46 n(H+)n(H 2 0) = [k 53 n(H 2 ) + (k 57 + k 5S )n(e)]n(H 2 + ) (A7) 

H 3 0+ : k 53 n(H 2 + )n(H 2 ) = (fc 59 + k 60 )n(e)n(H 3 O + ) (A8) 

: k 47 n{H + )n{0 2 ) = fc 6 in(e)n(0 2 f ). (A9) 
From (|A1|) - (|A9|) . the chemical equilibrium values are 

y(H-) = ^ (A10) 
k 3 



„, rm k 37 y(C)y(H) + k 28 y(C)y(H 2 ) + k 62 y(C)y(H 2 ) 

V{CH) = k 21 y(H) + h 2 y(0) (AU) 

y(CHA - k ^ C ^ H ^ ± k3iy(H 2 )y(CH) 
2 k 22 y{H) 

s(co+) . ^p^m 

k 54 y(H) 

(n+ , _ k 30 y(H+)y(O) + k 49 y(C+)y(Q 2 ) 
k 50 y(H) + k 5 iy(H 2 ) 

+ = hiy(Q + )y(H2) + k 4 ^y(H + )y(OH) 
Jy ' k 52 y{H 2 ) + k 56 y(e) V ' 

fnn+ s k 52 y(OH+)y(H 2 ) + k 46 y(H+)y(H 2 Q) 
k^ 3 y{H 2 ) + (fc 58 + k 57 )y{e) 

+ _ bgmgi (A17) 

(«60 + k 59 )y{e) 

Vi0t) - fc ^f + y) . (A18) 

By using the one-zone models, we have compared the results by the chemical model above and those by fully non- 
equilibrium calculations with the same set of reactions and have found that the evolution of temperature, as well as 
abundances of the coolants abundances agree excellently. 
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Fig. Bl. — The evolution of temperatures in prestellar cloud cores with metallicities Z/Zq = 0, 10 -6 , 10 -5 , 10 -4 , 10 -3 , 10~ 2 , 10 _1 , and 
1, as functions of the number density, which is calculated by one-zone models. The dashed lines indicate the constant Jeans masses. For 
those above 1O 2 M0 (below IMq), the gas is assumed to be fully atomic (molecular) in drawing those lines. 

B. CO AND OH COOLING FUNCTIONS 

We compute the cooling rate of CO and OH rotational transitions by counting the level populations under the 
assumption of statistical equilibrium (see Omukai 2001 Appendix B). The level energies and transitional coefficients 
are provided by Leiden Atomic and Molecular Database (Schoier et al. 2005). The dependence of the results for 
molecule M on the three variables, the H2 number density n(H2), column density parameter iV(M) and temperature 
T is fitted by the following four-parameter expression according to Neufeld & Kaufman (1993): 



1 

Lm 



1 



n(H 2 ) 
C 



LTE 



1 



n(H 2 ) 



H/2 



1 - 



-LTE 



(Bl) 



The column density parameter is defined as N(M) = N(M) / y ^(Jik) ( cm /k m s )j where iV(M) is the column 
density, and m(M) is the mass of molecule M. The fitting parameters Lo, which is a function of T alone, and £lte, 
ni/2, and a, which are functions of both N(M) and T, are presented in Table lB2l for CO and in Table lB3l for OH. 




Fig. B2. — The evolution of temperatures at the center of cloud cores during the prestellar collapse for various metallicities. This is 
calculated by one-zone model until 10 4 cm — 3 (dotted vertical line) and by hydrodynamical models for the higher density. The constant 
Jeans masses are indicated by the dashed lines. 
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FlG. B3. — The cooling and heating rates per unit mass by individual processes at the centers of the cloud cores as functions of number 
density. To avoid violent oscillatory behaviors and keep the compressional heating rate positive during the adiabatic phases, we plot them 
at every time when the density exceeds 1.1 times the previous values. 
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Fig. B4. — The evolution of the H2 concentrations y(ti2) at the centers of the cores. The state of y(H2) = 1/2 corresponds to the fully 
molecular gas. 
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Fig. B5. — The effective ratio of specific heat 7 = dlogp/dlogp at the center during the prestellar collapse for different metallicities as a 
function of number density. The horizontal lines show the critical values for fragmentation (1) and for hydrostatic core formation (4/3), 
respectively. 
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value for the isothermal Larson-Pcnston similarity solution (3.0) and the threshold for quasi-static collapse (1.0). 
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Fig. B7. — The radial distributions of (a) density (b) temperature (c) velocity and (d) enclosed mass at different epochs during the 
prestellar collapse for the cloud core with metallicity Z=0. The first curves shown correspond to the initial state of the calculation. They 
are plotted when the central density is enhanced by a factor of 10 2 ' 5 from the previous one except the last state shown, which corresponds 
to the time somewhat after the formation of a protostar. 
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Fig. B15. — The radial distributions of velocity at the protostar-formation epochs for the cloud cores with different metallicitics. 
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Fig. B16. — The mass accretion rates as functions of the protostellar mass derived by equation l |13j) . The horizontal line at dM*/dt = 
4 X 10~ 3 MQ/yr indicates the value above which the stationary accretion becomes impossible by the radiation force due to the accretion 
luminosity during the Kelvin-Helmholtz contraction. 




Fig. B17. — The evolution of protostellar radius during the accretion phase. The thick (dashed) lines represent the evolution with 
different metallicities of Z = 0, 1(T 5 Zq, 1CT 3 Zq, and ltT 1 Zq in descending order (Z = 1 Zq, 1(T 2 Zq, 10~ 4 Zq, and 1(T 6 Z Q in 
ascending order). The open circles on the lines indicate the epoch when the total energy production rate by deuterium burning reaches 
80% of the steady burning rate I/D,st = A^<5d; where 5jj is the energy available from the deuterium burning per unit gas mass. The filled 
circles mark the epoch when the total energy production rate by hydrogen burning exceeds 80% of the luminosity at the stellar surface L*. 
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F ig. B 18. — Radiative transfer effects on the central temperature evolution. The solid lines are the cases with radiativ e tr ansfer (same as 
Fig |B2l but enlarged for clarity), while dashed lines are the cases with continuum cooling rate reduced as in equation d 1 4D . Although the 
temperatures are slightly higher near optically thin to thick transition regime in the radiative transfer cases, their differences are small. 
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Fig. B19. — Radiative transfer effects on the envelope temperature. Th e solid lines are the cases with radiative transfer, while dashed 
lines are the cases with continuum cooling rate reduced as in equation Q14I . The temperatures are higher outside the first pr otoste llar cores 
for the radiative transfer cases. Note that these differences are more pronouced than tho se in the central temperatures (Fig |B18l) owing to 
heating from inner hot regions, which is not included in the simple method by equation l|14[l. 
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Fig. B20. — The growth and decline of elongation factors during the dust-cooling phase. The cases of metallicities [M/H]=-5, -4, -3, 
-2, and -1 are shown. For comparison, the cases [M/H]=-6, and -5 by the one-zone model is also presented. The dotted horizontal line 
indicates the critical value ~ 3 for fragmentation. 
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Fig. B21. — Heating effect on the ambient gas by radiation from the central protostar. The solid lines show the prestellar temperature 
T pre as a function of radius. The dashed lines show the dust temperature due to protostellar radiative heating 7h ca t = L/lGnar 2 , where L 
is the protostellar luminosity. The protostellar masses M* of those lines are indicated by numbers. The innermost point of a dashed line 
corresponds to the location of the rarefaction wave at 2M* , inside of which the matter falls at almost free-fall rate and the dynamics is not 
affected by protostellar heating. In cases with [M/H]< —3, the lines for 30Mq are not shown because the total core mass is less than 60Mq 
and the entire core has been swept by the rarefaction wave at M* = 30Mq. Similarly, the T^oat line at M* = 3Mq is not illustrated in 
the IZq case, where that at M, = IMq is shown instead. In the region where T pre is higher than Th ea t, the protostellar heating effect is 
negligible. Since the actual dust temperature T<j ust is determined by the balance between the cooling by its thermal radiation and heating 
not only by the protostellar radiation but also by the collisional coupling with gas, the dust temperature in the accretion phase roughly 
coincides with the higher of Theat or T prc . If T^ust > T pTe (i.e., T^ aat > T prc ), the gas temperature approaches T^ust (— ?heat m this case) 
by the collisional coupling. 



34 



K. Omukai et al. 



TABLE Bl 
Reduced Chemical Reactions 



Number 


Reaction 


1(H1) 


H+ + e -» H + 7 


2(H2) 


H + e — > H~ + 7 


3(H3) 


H-+H^H 2 + c 


4(H4) 


H 2 + H -> 3H 


5(H5) 


3H -> H 2 + H 


6(H6) 


2H + H 2 -)■ 2H 2 


7(H7) 


2H 2 2H + H 2 


8(H8) 


2H + grain -> H 2 


11(D1) 


D + H+ -> D+ +H 


12(D2) 


D+ + H-S.D + H+ 


13(D3) 


D + H 2 -> D + HD 


14(D4) 


D+ + H 2 ->■ H+ + HD 


15(D5) 


HD + H -> H 2 + D 


16(D6) 


HD + H+ H 2 + D+ 


21(M1) 


H + CH -> C + H 2 


22(M2) 


H + CH 2 -> CH + H 2 


25(M5) 


H + OH H 2 + O 


26(M6) 


H + H 2 -> OH + H 2 


27(M7) 


H + 2 -> OH + O 


28(M8) 


C + H 2 CH + H 


29(M9) 


O + H 2 -> OH + H 


30(M10) 


H+ + O — > 0+ + H 


31(M11) 


H 2 + CH -> CH 2 + H 


34(M14) 


H 2 + OH-> H 2 + H 


35(M15) 


20H -> H 2 + O 


36(M16) 


0H + C0^C0 2 + H 


37(M17) 


C + H -> CH + 7 


38(M18) 


C + OH^CO + H 


39(M19) 


c + o 2 -> CO + 


40(M20) 


O + H -> OH + 7 


41(Mzl ) 


20 — > O2 + 7 


a Q/"i\/roo~\ 

4Z(1V±ZZ ) 


/~i 1 v f V t I TT 

O + wri — ^ UU + n 




O 4- OH — i 4- FT 
w wrti r w 2 1 1 


lijl 1V1 ZiU 1 


H+ 4- OTT — s. OH+ 4- FT 


Aflf'tJIOpS 


tt+ -u TT^n tt.-,o+ 4- tt 

n "T HTU r LL2U n 


a '7/"^/^o'7^ 

4 / ( JVlZ ( J 


rl^ + U2 — ► U 2 + Ti 


/I Q /T\ /TOO \ 

4o(Mzo ) 


1 r^u \ nn+ 1 tj 
+ <Jri — ¥ L <J + ri 


a n /i\ /r on\ 


1 \ 0+ 1 (~^r\ 
Kj~ + U2 — > + 


OU^lvloU ) 


0+ _L FT _i. TT+ J- O 


r:i /l\/TQi \ 
01 l^lvlol ) 


0+ -1- TT ^ i. OTT+ -1- TT 

W ~r rl2 — r \Jci -\- n 


52(M32) 


OH+ +Ho -> HoO+ 4- H 

V_/ 1 1 \^ 112 ' Hz 1 -' |^ -l-L 


53(M33) 


H 2 0+ +H 2 -> H3O+ + H 


54(M34) 


CO++H->H++CO 


55(M35) 


C++c->C + 7 


56(M36) 


OH+ + e -»■ O + H 


57(M37) 


H 2 0+ +e -> OH + H 


58(M38) 


H 2 0+ +c -> O + H 2 


59(M39) 


H3O+ +e ->■ H 2 + H 


60(M40) 


H3O+ +e ->■ OH + 2H 


61(M41) 


Oj + c ->■ 20 


62(M42) 


H 2 + C ->■ CH 2 + 7 
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l°glO 


(AO 








tcm 


pcra 


turc 
















10K 


20K 


30K 


50K 


80K 


100K 


300K 


600K 


1000K 


1500K 


2000K 












old 


(Lo 


/ergs 


cm 3 














24.77 


24.38 


24.21 


24.03 


23.89 


23.82 


23.42 


23.13 


22.91 


22.63 


22.28 












l°Sio 


(^LTE 


/ PT"P"S 

/ 




b ) 








14 





21.08 


20.35 


19.94 


19.45 


19.01 


18.80 


17.81 


17.23 


16.86 


16.66 


16.55 


14 


5 


21.09 


20.35 


19.95 


19.45 


19.01 


18.80 


17.81 


17.23 


16.86 


16.66 


16.55 


15 


.0 


21.11 


20.37 


19.96 


19.46 


19.01 


18.80 


17.81 


17.23 


16.86 


16.66 


16.55 


15 


5 


21.18 


20.40 


19.98 


19.47 


19.02 


18.81 


17.82 


17.23 


16.87 


16.66 


16.55 


16 


.0 


21.37 


20.51 


20.05 


19.52 


19.05 


18.83 


17.82 


17.23 


16.87 


16.66 


16.55 


16 


.5 


21.67 


20.73 


20.23 


19.64 


19.13 


18.90 


17.85 


17.25 


16.88 


16.67 


16.56 


17 


.0 


22.04 


21.05 


20.52 


19.87 


19.32 


19.06 


17.92 


17.28 


16.90 


16.69 


16.58 


17 


.5 


22.44 


21.42 


20.86 


20.19 


19.60 


19.33 


18.08 


17.38 


16.97 


16.75 


16.63 


18 


.0 


22.87 


21.82 


21.24 


20.55 


19.95 


19.66 


18.34 


17.59 


17.15 


16.91 


16.78 


18 


.5 


23.30 


22.23 


21.65 


20.94 


20.32 


20.03 


18.67 


17.89 


17.48 


17.26 


17.12 


19 


.0 


23.76 


22.66 


22.06 


21.35 


20.71 


20.42 


19.03 


18.26 


17.93 


17.74 


17.61 












logi n 


(n-t /o 

V I/ 2, 


/cm - ' 1 ) 












14 


o 


3.29 


3.49 


3.67 


3.97 


4.30 


4.46 


5.17 


5.47 


5.53 


5.30 


4.70 


14 


5 


3.27 


3.48 


3.66 


3.96 


4.30 


4.45 


5.16 


5.47 


5.53 


5.30 


4.70 


15 





3.22 


3.45 


3.64 


3.94 


4.29 


4.45 


5.16 


5.47 


5.53 


5.30 


4.70 


15 


5 


3.07 


3.34 


3.56 


3.89 


4.26 


4.42 


5.15 


5.46 


5.52 


5.30 


4.70 


16 





2.72 


3.09 


3.35 


3.74 


4.16 


4.34 


5.13 


5.45 


5.51 


5.29 


4.68 


16 


.5 


2.24 


2.65 


2.95 


3.42 


3.92 


4.14 


5.06 


5.41 


5.48 


5.26 


4.64 


17 





1.74 


2.15 


2.47 


2.95 


3.49 


3.74 


4.86 


5.30 


5.39 


5.17 


4.53 


17 


.5 


1.24 


1.65 


1.97 


2.45 


3.00 


3.25 


4.47 


5.02 


5.16 


4.94 


4.27 


18 





0.742 


1.15 


1.47 


1.95 


2.50 


2.75 


3.98 


4.57 


4.73 


4.52 


3.84 


18 


5 


0.242 


0.652 


0.966 


1.45 


2.00 


2.25 


3.48 


4.07 


4.24 


4.03 


3.35 


19 


.0 


-0.258 


0.152 


0.466 


0.954 


1.50 


1.75 


2.98 


3.57 


3.74 


3.53 


2.85 


a 


14 


.0 


0.439 


0.409 


0.392 


0.370 


0.361 


0.357 


0.385 


0.437 


0.428 


0.354 


0.322 


14 


.5 


0.436 


0.407 


0.391 


0.368 


0.359 


0.356 


0.385 


0.437 


0.427 


0.354 


0.322 


15 





0.428 


0.401 


0.385 


0.364 


0.356 


0.352 


0.383 


0.436 


0.427 


0.352 


0.320 


15 


.5 


0.416 


0.388 


0.373 


0.353 


0.347 


0.345 


0.380 


0.434 


0.425 


0.349 


0.316 


16 


.0 


0.416 


0.378 


0.360 


0.338 


0.332 


0.330 


0.371 


0.429 


0.421 


0.341 


0.307 


16 


.5 


0.450 


0.396 


0.367 


0.334 


0.322 


0.317 


0.355 


0.419 


0.414 


0.329 


0.292 


17 


.0 


0.492 


0.435 


0.403 


0.362 


0.339 


0.329 


0.343 


0.406 


0.401 


0.317 


0.276 


17 


.5 


0.529 


0.473 


0.441 


0.404 


0.381 


0.370 


0.362 


0.410 


0.392 


0.316 


0.272 


18 





0.555 


0.503 


0.473 


0.440 


0.423 


0.414 


0.418 


0.446 


0.404 


0.335 


0.289 


18 


5 


0.582 


0.528 


0.499 


0.469 


0.457 


0.451 


0.470 


0.487 


0.432 


0.364 


0.310 


19. 





0.596 


0.546 


0.519 


0.492 


0.483 


0.479 


0.510 


0.516 


0.448 


0.372 


0.313 



TABLE B2 
Rotational cooling parameters for CO 
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logio(iV) 








temperature 








30K 


50K 


80K 


100K 


300K 


600K 






— loffi n 


(L 


/ergs cm a s~ 1 ) 








25.31 


24.53 


24.02 


23.82 


23.16 


22.90 






— l°gio 


(^LTE 


/ ergs cni^ s ^ ) 






10.0 


16.20 


15.44 


14.90 


14.67 


13.81 


13.58 


11.0 


16.20 


15.45 


14.90 


14.67 


13.81 


13.58 


12.0 


16.22 


15.46 


14.91 


14.67 


13.81 


13.59 


13.0 


16.36 


15.55 


14.96 


14.71 


13.83 


13.59 


14.0 


17.02 


16.00 


15.25 


14.95 


13.94 


13.67 


15.0 


17.74 


16.60 


15.77 


15.45 


14.49 


14.15 


16.0 


18.58 


17.38 


16.50 


16.16 


15.10 


14.77 


17.0 


19.42 


18.28 


17.40 


17.05 


15.84 


15.40 


18.0 


20.38 


19.21 


18.36 


18.02 


16.83 


16.35 






&1U 


(n-\ /o 

V 1/ z 


/cm - d ) 






10.0 


9.05 


8.98 


8.90 


8.88 


8.77 


8.72 


11.0 


9.05 


8.98 


8.90 


8.87 


8.77 


8.71 


12.0 


9.03 


8.97 


8.89 


8.87 


8.76 


8.71 


13.0 


8.88 


8.85 


8.81 


8.80 


8.74 


8.69 


14.0 


8.20 


8.27 


8.36 


8.40 


8.50 


8.50 


15.0 


7.25 


7.45 


7.71 


7.79 


7.99 


8.01 


16.0 


6.26 


6.46 


6.75 


6.85 


7.22 


7.29 


17.0 


5.26 


5.46 


5.75 


5.85 


6.23 


6.33 


18.0 


4.26 


4.46 


4.75 


4.85 


5.23 


5.33 


a 


10.0 


0.353 


0.543 


0.536 


0.530 


0.466 


0.434 


11.0 


0.354 


0.543 


0.536 


0.530 


0.466 


0.434 


12.0 


0.364 


0.544 


0.536 


0.529 


0.466 


0.434 


13.0 


0.399 


0.556 


0.534 


0.525 


0.463 


0.432 


14.0 


0.655 


0.586 


0.532 


0.521 


0.447 


0.424 


15.0 


0.641 


0.538 


0.538 


0.538 


0.448 


0.443 


16.0 


0.649 


0.572 


0.552 


0.516 


0.381 


0.370 


17.0 


0.566 


0.618 


0.580 


0.538 


0.381 


0.363 


18.0 


0.701 


0.645 


0.598 


0.550 


0.383 


0.369 



TABLE B3 
Rotational cooling parameters for OH 



